% PLOT_TREE   Plots a tree.
% (trees package)
%
% HP = plot_tree (intree, color, DD, ipart, res, options)
% -------------------------------------------------------
%
% plots a directed graph contained in intree. Many settings allow to play
% with the output results. Colour handling is differently on line plots than
% on patchy '-b' or '-p'. Even if metrics are nonexistent plot_tree will
% plot its best guess for a reasonable tree (see "xdend_tree"). Line plots
% are always slower than any patch display.
%
% Input
% -----
% - intree::integer:index of tree in trees or structured tree
% - color::RGB 3-tupel, vector or matrix: RGB values {DEFAULT [0 0 0]}
%     if vector then values are treated in colormap (must contain one value
%     per node then!).
%     if matrix (N x 3) then individual colors are mapped to each
%     element, works only on line-plots
% - DD:: 1x3 vector: coordinates offset {DEFAULT no offset [0,0,0]}
% - ipart::index:index to the subpart to be plotted {DEFAULT: all nodes}
% - res::integer>1: resolution for cylinders. Does not affect line and
%     quiver or blatt. {DEFAULT: 8}
% - options::string: {DEFAULT '-p'} has to be one of the following:
%     '-b'  : 2D pieces are displayed on a 3D grid (-b stands for -blatt)
%             showing the diameter but not as real cylinders. Output is a
%             series of patches. Fastest representation.
%        '-b1' : patches are mapped on x y {DEFAULT}
%        '-b2' : patches are mapped on x z
%        '-b3' : patches are mapped on y z
%     '-p'  : correct cylinder representation but not yet flawless and a
%             bit slower than "blatt" representation.
%     '-2l' : 2D (using only X and Y). forces line output (2D), no diameter
%             (slower), color is mapped independently of matlab, always
%             min to max.
%     '-3l' : 3D. forces line output (2D), no diameter (slower, as '-l')
%     '-2q' : 2D (using only X and Y). edges are represented as arrows
%             (using quiver) . Color vectors do not work.
%     '-3q' : 3D. edges are represented as arrows (using quiver, as '-q')
%   additional options:
%     '-thin'  : all diameters 1um, for line and quiver linewidth 0.25
%     '-thick' : all diameters + 3um, for line and quiver linewidth 3
%
% Output
% ------
% - HP::handles: depending on options HP links to the graphical objects.
%
% Example
% -------
% plot_tree (sample_tree)
%
% See also vtext_tree xplore_tree
% Uses cyl_tree ver_tree
%
% directly adapted for TREES toolbox, code from: Friedrich Forstner
%
% the TREES toolbox: edit, visualize and analyze neuronal trees
% Copyright (C) 2009  Hermann Cuntz

function HP = plot_tree (intree, color, DD, ipart, res, options)

% trees : contains the tree structures in the trees package
global trees

if (nargin < 1)||isempty(intree),
    intree = length(trees); % {DEFAULT tree: last tree in trees cell array}
end;

ver_tree (intree); % verify that input is a tree structure

% use full tree for this function
if ~isstruct (intree),
    tree = trees {intree};
else
    tree = intree;
end

if (~isfield(tree,'X'))||(~isfield(tree,'Y'))
    [xdend tree] = xdend_tree (intree); % if metrics are missing replace by equivalent tree
end

N = size (tree.X, 1); % number of nodes in tree

if (nargin < 4)||isempty(ipart),
    ipart = (1 : N)'; % {DEFAULT index: select all nodes/points}
end

if (nargin < 2)||isempty(color),
    color = [0 0 0]; % {DEFAULT color: black}
end;
if (size (color, 1) == N) && (size (ipart, 1) ~= N),
    color = color (ipart);
end

if (nargin < 3)||isempty(DD),
    DD = [0 0 0]; % {DEFAULT 3-tupel: no spatial displacement from the root}
end
if length(DD)<3,
    DD = [DD zeros(1, 3 - length (DD))]; % append 3-tupel with zeros
end

if (nargin < 5)||isempty(res),
    res = 8; % {DEFAULT: 8 points around cylinder}
end

if (nargin < 6)||isempty(options),
    options = '-p'; % {DEFAULT: full cylinder representation}
end

if strfind (options,'-b'),
    if isfield (tree, 'D'),
        D = tree.D (ipart);
    else
        D = ones (length (ipart), 1);
    end
    if strfind (options, '-thin'),
        D = ones (length (ipart), 1);
    end
    if strfind (options, '-thick'),
        D = D + 3;
    end
    [X1 X2 Y1 Y2 Z1 Z2] = cyl_tree (intree);
    X1 = X1 (ipart) + DD (1); X2 = X2 (ipart) + DD (1);
    Y1 = Y1 (ipart) + DD (2); Y2 = Y2 (ipart) + DD (2);
    Z1 = Z1 (ipart) + DD (3); Z2 = Z2 (ipart) + DD (3);
    warning ('off', 'MATLAB:divideByZero');
    % direction vectors:
    dP = [X2-X1 Y2-Y1 Z2-Z1] ./ repmat (sqrt ((X2 - X1).^2 + (Y2 - Y1).^2 + (Z2 - Z1).^2), 1, 3);
    warning ('on',  'MATLAB:divideByZero');
    a1 = strfind (options, '-b');
    if length (options) > a1 + 1
        typ = str2double (options (a1 + 2));
        if isnan (typ),
            typ = 1;
        end
    else
        typ = 1;
    end
    if isfield (tree, 'frustum') && (tree.frustum == 1),
        idpar = idpar_tree (tree);
        parD  = D (idpar);
    end
    switch typ % draw plates (4 coordinates back and forth on cylinder):
        case 2
            dP = dP (:, [1 3]);
            % use rotation matrix to rotate the data
            V1 = (dP*[0 -1;1  0]) .* (repmat (D, 1, 2) ./ 2);
            V2 = (dP*[0  1;-1 0]) .* (repmat (D, 1, 2) ./ 2);
            if isfield (tree, 'frustum') && (tree.frustum == 1),
                V3 = (dP*[0 -1;1  0]) .* (repmat (parD, 1, 2) ./ 2);
                V4 = (dP*[0  1;-1 0]) .* (repmat (parD, 1, 2) ./ 2);
                MX = [X1+V4(:, 1) X1+V3(:, 1) X2+V1(:, 1) X2+V2(:, 1)]';
                MY = [Y1          Y1          Y2          Y2]';
                MZ = [Z1+V4(:, 2) Z1+V3(:, 2) Z2+V1(:, 2) Z2+V2(:, 2)]';
            else
                MX = [X1+V2(:, 1) X1+V1(:, 1) X2+V1(:, 1) X2+V2(:, 1)]';
                MY = [Y1          Y1          Y2          Y2]';
                MZ = [Z1+V2(:, 2) Z1+V1(:, 2) Z2+V1(:, 2) Z2+V2(:, 2)]';
            end
        case 3
            dP = dP (:, [2 3]);
            % use rotation matrix to rotate the data
            V1 =(dP*[0 -1;1  0]) .* (repmat (D, 1, 2) ./ 2);
            V2 =(dP*[0  1;-1 0]) .* (repmat (D, 1, 2) ./ 2);
            if isfield (tree, 'frustum') && (tree.frustum == 1),
                V3 = (dP*[0 -1;1  0]) .* (repmat (parD, 1, 2) ./ 2);
                V4 = (dP*[0  1;-1 0]) .* (repmat (parD, 1, 2) ./ 2);
                MX = [X1          X1          X2          X2]';
                MY = [Y1+V4(:, 1) Y1+V3(:, 1) Y2+V1(:, 1) Y2+V2(:, 1)]';
                MZ = [Z1+V4(:, 2) Z1+V3(:, 2) Z2+V1(:, 2) Z2+V2(:, 2)]';
            else
                MX = [X1          X1          X2          X2]';
                MY = [Y1+V2(:, 1) Y1+V1(:, 1) Y2+V1(:, 1) Y2+V2(:, 1)]';
                MZ = [Z1+V2(:, 2) Z1+V1(:, 2) Z2+V1(:, 2) Z2+V2(:, 2)]';
            end
        otherwise
            dP = dP (:, [1 2]);
            % use rotation matrix to rotate the data
            V1 = (dP*[0 -1;1  0]) .* (repmat (D, 1, 2) ./ 2);
            V2 = (dP*[0  1;-1 0]) .* (repmat (D, 1, 2) ./ 2);
            if isfield (tree, 'frustum') && (tree.frustum == 1),
                V3 = (dP*[0 -1;1  0]) .* (repmat (parD, 1, 2) ./ 2);
                V4 = (dP*[0  1;-1 0]) .* (repmat (parD, 1, 2) ./ 2);
                MX = [X1+V4(:, 1) X1+V3(:, 1) X2+V1(:, 1) X2+V2(:, 1)]';
                MY = [Y1+V4(:, 2) Y1+V3(:, 2) Y2+V1(:, 2) Y2+V2(:, 2)]';
                MZ = [Z1          Z1          Z2          Z2]';
            else
                MX = [X1+V2(:, 1) X1+V1(:, 1) X2+V1(:, 1) X2+V2(:, 1)]';
                MY = [Y1+V2(:, 2) Y1+V1(:, 2) Y2+V1(:, 2) Y2+V2(:, 2)]';
                MZ = [Z1          Z1          Z2          Z2]';
            end
    end
    if (size (color, 1) > 1) && (size (color, 2) == 1),
        HP = patch (MX, MY, MZ, color');
    else
        HP = patch (MX, MY, MZ, 0);
        if (size (color, 2) == 3) && (size (color, 1) == 1),
            set (HP, 'facecolor', color); % map color
        end
    end
end

if ~isempty ([strfind(options, '-2') strfind(options, '-3')]),
    % if color values are mapped:
    if size (color, 1) > 1,
        if size (color, 2) ~= 3,
            if islogical (color),
                color = double (color);
            end
            crange = [min(color) max(color)];
            % scaling of the vector
            if diff (crange) == 0,
                color = ones (size (color, 1), 1);
            else
                color = floor ((color - crange (1)) ./ ((crange (2) - crange (1)) ./ 64));
                color (color <  1) =  1;
                color (color > 64) = 64;
            end
            map    = colormap;
            colors = map (color, :);
        end
    end
    if ~isempty ([strfind(options, '-2l') strfind(options, '-3l')]),
        if strfind (options, '-2l'),
            [X1 X2 Y1 Y2]       = cyl_tree (intree, '-2d');
            HP = line([X1(ipart) X2(ipart)]' + DD (1), [Y1(ipart) Y2(ipart)]' + DD (2));
        end
        if strfind (options, '-3l'),
            [X1 X2 Y1 Y2 Z1 Z2] = cyl_tree (intree);
            HP = line([X1(ipart) X2(ipart)]' + DD (1), [Y1(ipart) Y2(ipart)]' + DD (2),...
                [Z1(ipart) Z2(ipart)]' + DD (3));
        end
        if size (color, 1) > 1,
            for ward = 1 : length (ipart)
                set (HP (ward), 'color', colors (ward, :));
            end
        else
            set (HP, 'color', color);
        end
    end
    if strfind (options, '-2q'),
        [X1 X2 Y1 Y2]       = cyl_tree (intree, '-2d');
        HP = quiver (X1 (ipart) + DD (1), Y1 (ipart) + DD (2), ...
            X2 (ipart) - X1 (ipart), Y2 (ipart) - Y1 (ipart), 0);
        if size (color, 1) > 1,
            color = [0 0 0];
        end
        set (HP, 'color', color);
    end
    if strfind (options, '-3q'),
        [X1 X2 Y1 Y2 Z1 Z2] = cyl_tree (intree);
        HP = quiver3 (X1 (ipart) + DD (1), Y1 (ipart) + DD (2), Z1 (ipart) + DD (3),...
            X2 (ipart) - X1 (ipart), Y2 (ipart) - Y1 (ipart), Z2 (ipart) - Z1 (ipart), 0);
        if size (color, 1) > 1,
            color = [0 0 0];
        end
        set (HP, 'color', color);
    end
    if strfind (options, '-thin'),
        set (HP, 'linewidth', 0.25);
    end
    if strfind (options, '-thick'),
        set (HP, 'linewidth', 3);
    end
end

if strfind (options, '-p'),
    [X1 X2 Y1 Y2 Z1 Z2] = cyl_tree (intree);
    X1 = X1 (ipart) + DD (1); X2 = X2 (ipart) + DD (1);
    Y1 = Y1 (ipart) + DD (2); Y2 = Y2 (ipart) + DD (2);
    Z1 = Z1 (ipart) + DD (3); Z2 = Z2 (ipart) + DD (3);
    P  = [X1 Y1 Z1]; % location of starting nodes
    dP = [X2-X1 Y2-Y1 Z2-Z1]; % direction vector
    N  = size (X1, 1); % number of nodes in tree
    % initialize sampling angles, remove duplicate angles (0 and 2*pi)
    theta = linspace (0, 2 * pi, res + 1); theta = theta (1 : end - 1);
    % initialize storage arrays
    b1 = zeros (N, 3); b2 = zeros (N, 3);
    % calculate orthonormal vectors to the direction vector of each compartment
    for ward = 1 : N
        % singular value decomposition
        v = null (dP (ward, :)); %%% BOTTLENECK
        % orthogonal bases setup
        b1 (ward, :) = v (:, 2); b2 (ward, :) = v(:, 1);
    end;
    % replicate vectors and reshape
    b1 = repmat  (b1,  1, 2 * res);
    b1 = reshape (b1', 3, 2 * res * size (b1, 1))';
    b2 = repmat  (b2,  1, 2 * res);
    b2 = reshape (b2', 3, 2 * res * size (b2, 1))';
    
    % this is the first translation vector list to move the cylinder
    % terminals away from null
    dP = cat (2, zeros (N, 3 * res), repmat (dP (1 : end, :), 1, res));
    dP = reshape (dP', 3, N * 2 * res)';
    % second translation list to move each cylinder to it's final destination
    P  = repmat  (P (1 : end, :), 1, 2 * res);
    P  = reshape (P',             3, 2 * res * N)';
    
    % setup diameters:
    % use frustums for ratio (radius/parent) >$threshold and < 1/treshhold to
    % avoid strange visualization results of thin branches leaving
    % from bigger ones, else: use cylinders
    if isfield (tree, 'D'),
        D = tree.D;
    else
        D = ones (length (X1),    1);
    end
    if strfind (options, '-thin'),
        D = ones (length (ipart), 1);
    end
    if strfind (options, '-thick'),
        D = D + 3;
    end
    if isfield (tree, 'frustum') && (tree.frustum == 1),
        threshold = .15;
        idpar    = idpar_tree (intree); % vector containing index to direct parent
        parD     = D (idpar (1 : end));
        D        = D (ipart); parD = parD (ipart);
        use_parD = double ((D (1 : end) ./ parD) > threshold &...
            (D (1 : end) ./ parD) < (1 / threshold));
        diff_rad = parD - D (1 : end);
        D        = cat (2, repmat (D (1 : end) + (use_parD .* diff_rad), 1, res),...
            repmat (D (1 : end), 1, res));
        D        = reshape (D',  1, res * 2 * N)';
    else
        D        = D (ipart);
        D        = repmat  (D,   1, res * 2); D = reshape (D', numel (D), 1);
    end
    % replicate angular sampling
    theta_sin = repmat (sin (theta)', N * 2, 3);
    theta_cos = repmat (cos (theta)', N * 2, 3);
    % create final vertex list
    vertex_array = dP + P + (theta_cos.*b1 + theta_sin.*b2) .* repmat (D/2, 1, 3);
    % face index creation, start with the two points on the base level and
    % then connect them to the two others on the terminal level
    % start with a standard poly list for a cylinder then replicate and shift
    first_col  = (1       : 1 : res - 1)';
    second_col = (2       : 1 : res)';
    third_col  = (res + 2 : 1 : 2 * res)';
    fourth_col = (res + 1 : 1 : 2 * res - 1)';
    poly_array = [first_col, second_col, third_col, fourth_col];
    poly_array = [poly_array; [res, 1, res+1, 2*res]];
    poly_array = repmat  (poly_array, N, 1);
    shift_vec  = repmat  ((0 : 2 * res : (2 * res * N) - 1)', 1, res * 4);
    shift_vec  = reshape (shift_vec', 4, res * N)';
    poly_array = poly_array + shift_vec;
    if (size (color, 1) > 1) && (size (color, 2) == 1),
        C = repmat (color, 1, res * 2); C = reshape (C', numel (C), 1);
        HP = patch ('Faces', poly_array, 'Vertices', vertex_array, ...
            'cdata',          C', 'facecolor', 'interp', 'linestyle', 'none');
    else
        C = repmat (color, size (vertex_array, 1), 1);
        HP = patch ('Faces', poly_array, 'Vertices', vertex_array, ...
            'faceVertexCData', C, 'faceColor', 'interp', 'linestyle', 'none');
    end
end

if ~(sum (get (gca, 'DataAspectRatio') == [1 1 1]) == 3)
    axis equal
end
